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COMPUTER PROGRAMS FOR THERMODYNAMIC 
AND TRANSPORT PROPERTIES OF HYDROGEN 
(Tabcode-II) 


H. M. Rodar, R. D. McCarty, and W. J. Hall 

The thermodynamic and transport properties of para and 
equilibrium hydrogen have been programmed into a series of com- 
puter routines. Input variables are the pair's pressure- temperature 
and pressure-enthalpy. The programs cover the range from 1 to 
5000 psia (34 MN/m^) with temperatures from the triple point to 
6000° R (3300 K) or enthalpies from -130 BTU/lb (-623 J/mol) to 
25,000 BTU/lb (117000 J/mol). Output variables are enthalpy or 
temperature, density, entropy, thermal conductivity, viscosity, 
velocity of sound, heat capacity at constant pressure, heat capacity 
at constant volume, the heat capacity ratio, and a heat transfer 
parameter. Property values on the liquid and vapor boundaries are 
conveniently obtained through two small routines. The programs 
achieve high speed by using linear interpolation in a grid of precom- 
puted points which define the surface of the property returned. The 
maximum errors arising from the linear interpolation are shown on 
individual deviation plots for each combination of variables. Error 
estimates for the sources of data are similarly displayed. 

Key Words: Computer programs; density; enthalpy; entropy; 
heat capacity at constant pressure; heat capacity at constant 
volume; heat capacity ratio; heat transfer coefficient; hydrogen; 
pressure; saturation boundary; temperature; thermal conductivity; 
velocity of sound; viscosity. 

1. INTRODUCTION 

Extensive use of hydrogen in the chemical industry and as a propellant in the U. S. 
space program requires that the engineer be provided values of thermodynamic and transport 
properties of this fluid in useful form. In many problems the properties have to be evaluated 
at numerous thermodynamic state points. The problems are often handled on a computer 
and they involve iterative computational methods. In the present report we describe a 
series of computer routines which are fast, moderately accurate, and operate over wide 
ranges of input variables. They are based on selected source data. The input (entry) and 
output (returned) variables are: 

In P ut Oitput Input Output 

pressure- enthalpy pressure- temperature, and phase or state 

temperature density enthalpy density 

entropy entropy 

thermal conductivity thermal conductivity 

viscosity viscosity 

velocity of sound velocity of sound 

heat capacity, C heat capacity, C 

heat capacity, heat capacity, 

heat capacity ratio heat capacity ratio 

heat transfer parameter heat transfer parameter 









i 


\ 


In addition there are two small saturation routines which return values on the satura- 
tion boundaries. The computer routines were developed primarily for nuclear rocket engine 
design. The programs cover a wide range of conditions as follows: for temperatures 

from the triple point (24.845°R, 13.803 K) to about 6000°R (3300 K), for pressures from 3 to 
2 

5000 psia (34 MN/m ), and for enthalpies from - 130 to 25,000 BTU/lb (-623 to 117000 J/mol). 
In addition, the thermodynamic and transport properties of both para hydrogen and equilibrium 
hydrogen are available through these programs. An important criterion for many applications 
is computational speed, as some process control calculations have to be made in real time. 
The linear interpolation employed achieves speed but sacrifices accuracy. As an aid to the 
user, we have determined the maximum error due to the interpolation scheme and estimated 
the error in the source data. These errors are displayed in departure graphs for each 
program. 

An older version of the para hydrogen programs, which has been distributed widely, 
should now be considered obsolete (Tabcode-I, Hall, et al. , 1967, and McCarty 1968). The 
changes and improvements over the older version are as follows: 

1. The temperature range has been extended to 6000° R, the enthalpy range to 25,000 
BTU/lb. 

2. New measurements of thermal conductivity between 30° and 275° R are now included. 

3. Improved values for both thermal conductivity and viscosity have been inserted for 
temperatures above 1200 P R. 

4. The saturation subroutines are now included in the package. 

5. For equilibrium hydrogen, which was not treated previously, an entire set of programs 
is now provided. 

6. A combination of variables is now included which describe a commonly used heat 
transfer parameter, designated herein as L-factor (see Section 1. 1). 

The main divisions of the report are text, program listings, and error graphs. The 
text contains all the normal description given for computer programs such as deck size, 
calling sequence, etc. However, for this problem the selection of source data and the com- 
putation of the property values placed in data arrays are important; they are, therefore, 
described in some detail. 


1. 1 Symbols and Units Used 

The programs were developed for engineering applications; therefore, the units of 
this report are those commonly used in engineering, i.e., psia, BTU, etc. The sources 
of data, however, normally employ some version of metric units. In the discussion of the 
equation of state, the virial B, and in figures 11 and 12, the units of the original sources 
had to be retained. Conversion factors for changes from source papers to generating 
functions to the present program have all been taken from NASA Special Publication 7012 
(Mechtly 1969). 
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pressure, psia 
temperature, °R 
enthalpy, BTU/lb 

* a enthalpy of saturated vapor, BTU/lb 
ga s 

^liquid enthal Py of saturated liquid, BTU/lb 
entropy, BTU/lb - °R 
density, lb/ft 3 

specific heat at constant pressure, BTU/lb - °R 

specific heat at constant volume, BTU/lb - °R 

the ratio C /C , dimensionless 
P v 

sound velocity, ft/s 

viscosity, lb-h/ft a 

thermal conductivity BTU/h-ft - °R 

a heat transfer parameter equal to k 0 - 6 C 0 , 4 /t) 0 - 4 , with units as 

P 

shown above. [ This heat transfer parameter is the combination of 
thermal properties occurring in the conventional Dittus-Boelter 
(Nusselt type) equation describing the forced convection heat 
transfer coefficient for turbulent flow (see, for example, 
equation 9- 10a, page 219, McAdams 1954)]. 
the phase designator, dimensionless 


Additional symbols used in text : 

B or B(T) the second virial coefficient 


superscript 


ideal gas internal energy at absolute zero, the reference state, 

taken to be zero 

general symbol for property 

critical pressure 

gas constant 

general symbol for the return from the liquid-vapor two phase region 
ideal gas thermodynamic properties, i. e. , C°, C°, H° 

p’ v* 

internal energy 

volume 

PV/RT 


Important state pointa for para and equilibrium hydrogen; 


triple point 
P = 1. 022 naia 
T = 24. 845*R 


normal boiling point 
P= 14.696 psia 
T = 36. 482* R 


critical point 
P = 187. 506 psia 
T = 59. 356*R 
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2. HYDROGEN MODIFICATIONS, PROPERTY DIFFERENCES, AND PHASE DIAGRAMS 

A brief review of the more important features in the physics of hydrogen is essential 
if we are to understand what properties the programs return, and why some of the 
peculiarities of the programs occur. 

2. 1 Hydrogen Modifications 

Hydrogen molecules occur in two natural modifications, parahydrogen and ortho- 
hydrogen. In the parahydrogen molecule the two nuclear spins are opposed while in ortho- 
hydrogen they are alined. Under normal conditions, the transition probability between 
states is practically zero; a catalyst, a sharp gradient in magnetic field, or external radia- 
tion (e.g. , nuclear radiation, x-rays, ultraviolet, etc.) is required to let the transition pro- 
ceed. Conversion between these two modifications is temperature dependent, as shown in 
figure 1, which is based on Woolley, et al. , (1948). At low temperatures in the oresence of 
an appropriate catalyst, the thermal equi* >rium favors the para modification. At the 
normal boiling point (36.482°R), the equ .oriurn composition is 99.79% para and only 0.21% 
orthohydrogen. This composition is often simply referred to as parahydrogen or p-hydrogen. 


TEMPERATURE, K 



Figure 1. Hydrogen Composition at Equilibrium 
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Conversely, the equilibrium composition near ambient temperature is 75% ortho and 25% 
para. This latter composition is usually designated normal or n-hydrogen. When conversion 
occurs between the normal boiling point and room temperature, then the thermodynamic 
equilibrium mixture of para and ortho hydrogen is called "equilibrium" or e-hydrogen and 
the composition at intermediate temperatures varies as shown in figure 1. In this report we 
describe programs for both p-hydrogen and e-hydrogen. 

At temperatures above 2700°R (1500 K) the hydrogen molecule dissociates into hydro- 
gen atoms. This dissociation is both temperature and pressure dependent. Assumptions 
made in this report are that dissociation occurs according to the equilibrium constant for 
normal hydrogen, and that the property values can be obtained by adding the dissociation 
effects calculated for the ideal gas to the extrapolated properties of the real gas. 

2.2 Property Differences 

At comparable state points, property values of the various modifications of hydrogen 
may differ, that is, some of the properties depend on ortho-para composition. Experi- 
mentally these effects are observed in different vapor pressures and in different critical 
parameters. Most important to this discussion is the fact that the properties of the ideal 
gas for para and equilibrium hydrogen as calculated from statistical mechanics differ, and 
further, since the composition change depends only on temperature, that the property 
differences depend only on temperature. Thus, if a difference exists for an ideal gas prop- 
erty at low pressures, then the same difference exists at all other pressure levels for a 
given temperature. Stated another way, the statistical calculation, by including the temper- 
ature dependent composition change, automatically includes the reaction occurring as para- 
transforms to equilibrium-hydrogen — i. e. , the sensible heat of e-hydrogen includes the 
appropriate part of the heat of chemical reaction. Not all properties are cor. • osition depen- 
dent; density and viscosity, for example, are nearly independent of composition. On the 
other hand, the properties with significant ortho-para dependency are specific heat, and 
properties related to specific heat such as velocity of sound, entropy, enthalpy, and thermal 
conductivity. When the properties of p- and e- hydrogen differ, they may differ considerably. 
For example, the value of enthalpy near the critical point for p-hydrogen is about 16 BTU/lb, 
whereas the value for equilibrium hydrogen is about 31 BTU/lb. The difference in prop- 
erties for the ideal gas specific heats, C • , is illustrated in figure 2. Differences for the 
other properties are best expressed in terms of a table; approximate differences can be 
obtained by comparing corresponding entries for a pressure of 1 psia in Appendix A (pages 
38 and 45). 
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figure 2. Ideal Gas Specific Heat of Para and I^quilibrium Hydrogen 


2. 3 Range of the Programs and Phase Diagrams 

For programs which use pressure and temperature as entry variables, the pressure 
range is 1 to 5000 psia and the temperature range is from the triple point (24.845°R) to 
6000° R. The P-T range of the programs for both para and equilibrium hydrogen is illus- 
trated by the shaded area in figure 3. The programs cover property values tor liquid, gas or 
vapor states, and dissociation. Note that the programs do not cover the solid phase nor 
values at pressures less than 1 psia. The low pressure limit is arbitrary; it arises from the 
behavior of entropy as follows. At any given temperature, the value of entropy approaches 
infinity as the pressure approaches zero. Since we cannot represent an infinity with linear 
interpolation between fixed points, a cutoff, 1 psia, was chosen. 

For the pro t s which use pressure and enthalpy as entry variables, the pressure 
range is the same, 1 i ia to 5000 psia. The range of enthalpy is -130 BTU/lb to 25,000 
BTU/lb. The enthalpy limit3 correspond to the enthalpy of the liquid at the triple point, and 
the ideal gas enthalpy at the highest temperature, 6000 R. The phase diagram of para- 
hydrogen in pressure- enthalpy coordinates is shown in figure 4. The phase diagram for 
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Figure 4. The P-H Phase Diagram for Parahydrogen 
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e-hydrogen would look quite similar to figure 4 except that the liquid- vapor two phase bound- 
ary is gradually displaced to higher values of enthalpy. For example, the enthalpy of 
e-hydrogen at the critical pressure is about 31 BTU/lb while the maximum value of enthalpy 
along the two-phase boundary is about 110 BTU/lb. The programs cover property values for 
liquid, gas or vapor states, and dissociation as before. They do not cover the solid phase, 
pressures below 1 psia, or the two phase liquid- solid ("slush") conditions. 

A comparison of figures 3 and 4 discloses that the conditions, the phases or states 
of the hydrogen, covered by the P- T and the P-H programs are similar, except that the P-H 
programs cover the range of states where liquid and gas coexist while the P-T programs do 
not. A concise designation of the phase or state of any given point is therefore available 
only by reference to P-H coordinates. Since the description of phase is important in many 
problems, a scheme to designate the phase was devised and is presented in the next section. 


2.4 The Phase Designator, Q 

The various physical states possible in the P-H phase diagram are described by the 
phase designator, Q, available only from program PHTEMP (see also page 11). Values of 
Q are shown in figure 4 and in table 1. Note that figure 4 includes scale changes at 
100 BTU/lb and at 200 psia. For the special case of the two phase liquid-vapor region, Q 
corresponds to the conventional quality, as defined on page 11. 

Table 1. Values of the Phase Designator, Q 


Value of Q 

Applicable portions of :he phase diagram, P = 187. 506 psia 

c 

1.0 

single phase gas or vapor 

single phase compressed fluid at pressures above P 
saturated vapor line C 

0.0 < Q < 1.0 

two phase liquid-vapor region 

0.0 

saturated liquid line 

compressed liquid (fluid) states at pressures below P 
liquid- solid boundary line C 

-1.0 < Q < 0.0 

two phase liquid- solid region 

-1.0 

solid - compressed solid boundary 

Q < -1.0 

compressed solid states 
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3. OVERVIEW OF THE PROGRAMS 
The decks available and described in this report are listed in table 2. 


Table 2. List of Available Program Decks 


Input 

Output 

p-hydrogen 

e-hydrogen 

P-T 

enthalpy- 

PTENTH 

ETENTH 


entropy 

TPENTR 

ETENTR 


density 

PTDENS 

ETDENS 


Cp 

PTCP 

ET CP 


Cv 

entry PTCV 

entry ETCV 


gamma 

entry PTGAMM 

entry ETGAMM 


sound velocity 

PTSOUN 

ETSOUN 


viscosity 

i'TVISC 

ETVISC 


thermal conductivity 

PTCOND 

ETCOND 


heat transfer parameter 

PTLFAC 

ETLFAC 

P-H 

temperature , and Q 

PHTEMP 

EHTEMP 


entropy 

PHENTR 

EHENTR 


density 

PHDENS 

EHDENS 


S 

PHCP 

EHCP 


Cv 

entry PHCV 

entry EHCV 


gamma 

entry PHGAMM 

entry EHGAMM 


sound velocity 

PHSOUN 

FHSOUN 


viscosity 

PH VISC 

EHVISC 


thermal conductivity 

PHCOND 

EHCOND 


heat transfer parameter 

PHLFAC 

EHLFAC 

Saturation Routines 



P vap 1 

T vap' H gas> H liquid 

PSATH 

EPSATH 

T vap j 

P H H 

vap' gas' “liquid 

TSATH 

ETSATH 


Including the small saturation routine, there are 36 decks (the output options C , C , 
and the ratio C p /C y are combined into one deck). Each deck exists in a CDC and an IBM 
version. Program listings are given in Appendix B. 
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3 . 1 Method U sed 
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The conceptual idea of the computer programs follows from the use of a T-S or 
Molliei r iag ram. Using such a chart to obtain property data one visually or graphically 
interpolates between constant property lines. The computer programs imitate this procedure 
on the computer. However, in contrast to th^ visual procedure where any combination of 
existing variables can be used as entries to the chart, the computer programs are restricted 
to tw^ entry possibilities, namely pres sure- temperature and pr e ssure- enthalpy . 

Each property forms a surface with the two entry variables. A pre-deter mined 
number of points on this surface are calculated and form the array of data. The surface is 
divided into regions, the boundaries of which depend upon the relative linearity of the partic- 
ular property surface. The P, T, or H spacings vary from region to region to compensate 
for the non- linear ity of the property surface. In general, spacings will be small where the 
property is varying quite rapidly and they will be large where the property is varying more 
slowly. Each program consists of an array of data and an interpolation package. The 
interpolation package finds the four points which bracket the entry values and then performs 
a two dimensional linear interpolation to obtain the value of the property to be returned. 

The arrays of data are necessarily limited. When coupled with linear interpolation, 
errors in the values returned are possible. We have, therefore, determined the error 
envelope for each program and these are displayed as departure graphs. 

3.2 Peculiarities of the P-T Programs 

One prominent feature shown in the appropriate phase diagram, figure 3, is the vapor 
pressure curve. In all of the P-T programs this line is mathematically thin, that is, the 
property value returned corresponds to either the liquid or the vapor side. In following an 
isobar, an isotherm, or any other path that crosses the vapor pressure curve, one will 
encounter sizeable changes in the different variables. In addition, for many variables these 
large changes also exist for conditions near or just above the critical point. Re- stating the 
peculiarity a different way: it is impossible to pick a P-T input which will place the coordi- 

nates inside a two phase region. Since it is often necessary to obtain property values for 
conditions of saturation, saturated liquid, saturated vapor, or both, we have provided a way 
of obtaining these values. The saturation routines listed in table 2 are designed to fill this 
need. They are described in Section 3.4. 

3.3 Peculiarities of the P-H Programs 

In contrast to the P-T programs it is quite possible to choose input values of P and H 
which place a point in question into a two phase region. From figure 4 it is seen that the 
boundaries of the programs are set specifically to omit returns from the liquid- solid (slush) 
two phase region. Returns from the liquid-vapor two phase region are, however, possible. 
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Three distinct possibilities or problems connected with these returns are discussed in seme 
detail below. 

The Phase Designator Q, from PHTEMP - One problem often encountered is that we need 
to know where we are on the phase diagram. To solve this problem the program PHTEMP 
(EHTEMP) is structured slightly different than all of the other programs. The other P-H 
programs carry the two entry variables P and H in the call list. PHTEMP carries an 
additional variable, Q, in its call list. A call to PHTEMP (P, H, Q) and a check of the Q 
that is returned against table 1 will locate the point on the phase diagram. For example, 
if 0 < Q < 1. then the P-H input is in the liquid-vapor coexistence region. For this region 
Q corresponds to the conventional quality and is defined as the percent (by mass) of the 
vapor phase present, or 


Q 


(H ’ H v -J / 
liquid 


H ’ H.. .J 
gas liquid 


where H^^ and H gas are the enthalpy values at the respective phase boundary appropriate 
to the entry pressure (see also fig. 4). As noted before, values of Q for other parts of the 
phase diagram are given in table 1. 

Ret urns from the Liquid- Vapor Two Phase Region - Assume that the input values of P and H 
place the point in question into the liquid-vapor two phase region. For this region the entry 
enthalpy H is used in all PH (EH) programs to find the phaue designator or quality as 


defined above, and the phase designator or quality is then used to define the return, RE. 
For temperature, density, and entropy, these returns have the conventional meaning which 
is given in table 3. Other variables, that is k, T), a>, C , and L, are not defined for this 
region in the sense of thermodynam.es. However, the engineer often uses these properties 
in this region, and then defines them according to the specific application. Therefore, the 
returns given are artificial involving Q and property values corresponding to the saturated 
liquid and saturated vapor as defined in table 3. C p and y are infinite in this two phase 
region; thus the largest value available to the machine is returned. The returns are 
summarized in table 3. 


Property Values on the Saturation Boundaries - Property values of saturation are required in 
many applications. From the preceding discussions it is clear that these values are simply 
a special case of returns from the liquid-vapor two phase region. A property value on a 
phase boundary is obtained from the appropriate PH (EH) program by using pressure and the 
saturation enthalpy, liquid or gas, as entries. The saturation enthalpies can be obtained 
from the two small saturation routines. 

3.4 The Saturation Routines 


For points on the saturation boundary, i.e. f on the vapor pressure curve of figure 3, 
either P or T is sufficient to specify the conditions of the point in question (phase rule). 
Obviously, the pressures have to be between 1.022 and 187.506 psia and the temperatures 


t 


0 



11 


- 4 . 


V 



Table 3. Returns from the Liquid- Vapor Two Phase Region 


Program 

return mode 

description of return 

PHTEMP 

EHTEMP 

conventional 

RE = the temperature corresponding 
to the input pressure (vapor pressure). 
Additional output: the phase designator 
0 . 

PHDENS 

EHDENS 

conventional 

RE = average density of mixture 

= 1 . 0/[Q/p + (1 -Q)/ p .,] 

L vapor liquid J 

PHENTR 

EHENTR 

conventional 

RE proportional to quantities of 

PHCOND 

EHCOND 

artificial 

liquid and vapor present according 

PHVISC 

EHVISC 

a r tificial 

to the equation 

PHSOUN 

EHSOUN 

artificial 


PHCV 

EHC V 

artificial 

RE - Q(- - tt )rn 

vapor liquid liquid 

PHLFAC 

EHLFAC 

artificial 

where v stands for the property 
in question 

PHCP 

EHCP 

infinite 

RE - the largest value available to 

PHGAMM 

EHGAMM 

infinite 

the machine 


between 24.845 and 59t356°R, The saturation routines need tr be entered ' ith only one 

variable; the corresponding other variable on the vapor pressure curve is returned. In 

addition, the enthalpies of the saturated liquid, H, . . and the saturated vapor, H , are 

liquid r gas 

returned corresponding to the two phase boundaries of figure 4. Saturation values for all 
other properties can then be obtained by using pressure and the appropriate saturation 
enthalpy as input to the desired P-H program, as shown in the example on page 16. 

4. PROGRAMMING INFORMATION 
4. 1 Deck Description 

All of the decks are coded as FORTRAN functions. The decks are completely 
self-contained, including the data arrays which are used for interpolation. Each deck may 
be used separately or in conjunction with others. 

A prefix PT (ET) on a subprogram name means that it requires pressure-temperature 
as input variables while the prefix PH (EH) indicates that the input variables are pressure- 
enthalpy. The remainder of the deck name designates the variable returned. The names of 
the decks and units for input and output are given in table 4 for p- hydrogen and in table 5 for 
e -hydrogen. Also shown in these tables are data loading routines, storage requirements, 
and number of cards in each deck. 
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Table 4. Deck Description, p-hydrogen 


Input variables 
and units 

Deck 

# 

Sub-program 
na me 

Output variable 
and unit 

Extra sub- 
routine for 

Storage 

requirement 

Number 
c; v_a rds 





data loading 


in deck 

pressure psia 

1 

PTENTH 

enthalpy 

cDC - 2 

1717 

CDC 25 8 

and 

temperature 



BTC/ Id 



IBM 25 1, 

* R 

2 

PTDENS 

density 

CD C - 3 

1460 

CDC 234 




lb/ ft 5 



IBM 227 


3 

PTENTR 

entropy 

CDC - 1 

1274 

C DC 208 




BTC/lb R 



IBM 2 08 


4 

PTCOND 

the rmal 

CLK. - i 

t i = 





conductivity 
BTC ft/ h ft 2 ' R 



IBM i -2 


5 

PTVISC 

vi scosity 

CDC - l 

978 

cdc i r 




lb h /ft* 



IBM l 9-. 


6 

PTSOUN 

velocity ot 

CDC - 1 

980 

CDC 164 




sound ft/s 



IBM 105 


7 

PTC P 

C , BTU/lb® R 
C p , BTU/lb® R 

CDC - 4 

2151 

CDC 20 2 



PTC V 

IBM - 1 


IBM 291 



PTGAMM 

V * NONE 




8 

PTLFAC 

P 

CDC - 1 

1 157 

CDC 18 3 
IBM 180 

pressure psia 
and 

3 

PHTEMP 

temperature 

°R 

quality see text 

CDC - 1 

917 

CDC 160 

enthalpy 

BTU/lb 





IBM 161 


10 

PHDENS 

density 

CDC - l 

1142 

CDC 193 




lb/ ft 3 



IBM 194 


11 

PHENTR 

entropy 

CDC - 1 

1130 

CDC 181 




B TU/lb" R 



IBM 180 


12 

PHCOND 

thermal 

CDC - 1 

1304 

CDC 213 




conductivity 
BTU ft/h ft 2 1 R 



IBM 209 


13 

PHVISC 

viscosity 

CDC - l 

1402 

CDC 266 




lb h /ft 2 



IBM 263 


14 

PHSOUN 

velocity of 

NONE 

795 

CDC 145 




sound, ft/s 



IBM 148 


15 

PHCP 

C , BTU/lb’R 
C P , BT(J/ lb* R 

CDC - 4 

2077 

CDC 294 



PHC V 

IBM - 1 


IBM 289 



PHGAMM 

Y, V NONE 




16 

PHLFAC 

k O.« C 0.4 /f1 0.4 

P 

CDC - 1 

1 304 

CDC 215 
IBM 210 

temperature 

17 

TSATH 

vapor pressure 

NONE 

161 

CDC 29 

*R 



H gas* BTU/lb 
H Uq( BTU/lb 



IBM 32 

pressure 

psU 

18 

PSATH 

sat. temperature 

NONE 

161 

CDC 29 



H gas* BTU/lb 
H Ur BTU/lb 



IBM 32 




f 


\ 
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Table 5, Deck Description, e-hydrogen 


Input variables 
and units 

Deck 

# 

Sub-program 

name 

Output variable 
and unit 

Extra sub- 
routine for 
data loading 

Storage 

requirement 

Number 
of cards 
in deck 

pressure psia 
and 

1 

ETENTH 

enthalpy 

BTU/lb 

CDC - 1 

1717 

CDC 279 
IBM 255 

temperature 
° R 


2 

ETDENS 

density 
lb /ft 3 

CDC - 2 

1460 

CDC 229 
IBM 225 


3 

ETENTR 

entropy 
BTU/lb R 

CDC - 1 

1274 

CDC 209 
IBM 206 


4 

ETCOND 

thermal 

conductivity 
BIT ft/ h ft 2 R 

CDC - I 

1 1 56 

CDC 184 
IBM 182 


5 

ETVISC 

viscos lty 
lb h /ft* 

CDC - 1 

978 

CDC 197 
IBM 195 


6 

ETSOUN 

velocity of 
sound ft/ s 

CDC - 1 

980 

CDC 164 
IBM 162 


7 

ETC P 
ETCV 
ETGAMM 

C p , BTU/lb R 
C v , BTU/lb R 
v, NONE 

CDC - 3 
IBM - 1 

2151 

CDC 291 
IBM 293 


8 

ETLFAC 

k°* P C** 4 /— 4 
P 

CDC - 1 

1157 

CDC 151 

IBM 180 

pressure psia 
and 
enthalpy 

9 

EHTEMP 

temperature 

J R 

quality see text 

CDC - 1 

917 

CDC 161 
IBM 160 

BTU/lb 


10 

EHDENS 

density 
lb /ft' 1 

CDC - 1 

1 142 

CDC 195 
IBM 194 


l 1 

EHENTR 

entropy 
BTU/lb 0 R 

CDC - \ 

1130 

CDC 182 
IBM 180 


12 

EHCOND 

thermal 
conductivity 
BTU ft / h ft 2 ° R 

CDC - 1 

1304 

CDC 213 
IBM 211 


13 

EHVISC 

viscosity 
lb h /ft 2 

CDC - 1 

1402 

CDC 267 
IBM 267 


14 

EHSOUN 

velocity of 
sound, ft/s 

NONE 

7°8 

CDC 141 
IBM 141 


15 

EHCP 

EHCV 

EHGAMM 

C BTU/ lb° R 
Cy, BTU/lb°R 
V, NONE 

CDC - 2 
IBM - l 

2077 

CDC 289 
IBM 291 


16 

EHLFAC 

c o, </T1 n.4 

P 

CDC - l 

1304 

CDC 227 
IBM 210 

temperature 

•R 

17 

ETSATH 

vapor pressure 
H BTU/lb 

H hq , BTU/lb 

NONE 

161 

CDC 30 
IBM 30 

psia 

18 

EPSATH 

sat, temperature 
H BTU/lb 

Hf iq . BTU/lb 

NONE 

161 

CDC 30 
IBM 30 
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Since nearly every computer has available some type of FORTRAN processing 
language, the programs were written in FORTRAN IV. Program listings for both the CDC 
and the IBM versions of all decks are given in Appendix B. 

4. 2 Calling Sequence 

Since the decks are individual FORTRAN functions the calling sequences are 
predictably simple. 

P = x. x 

T = y. y 

H = z. z 

and then 

RHO = PTDENS (P, T), 
or RHOE = EHDENS (P, H). 

Provided that the values of P, T. and H are valid, that is, within the program range, 

the density of para hydrogen is stored into RHO or the density of equilibrium hydrogen 
into RHOE. 

The only variation on this theme occurs for the routine PHTEMP (EHTEMP) and 
the saturation routines. The change occurs not in the calling sequence but rather in what 
information is available after the call is made. Examples are given below. 

P = X. X 
H = z. z 

TEMP = PHTE.ViPfP, H, Q) or TEMP = EHTEMPfP, H, Q). 

The temperature corresponding to the input pressure is found in TEMP. In addition, 
the phase designator Q is available through the parameter list, for example, 

QUAL = Q. 

Note that we did not set Q prior to the call. The saturation routines use very similar logic. 
P = X. X 

TVAP r PSATH(P, HG, HL) or TVAP = EPSATH(P, HG, HL) 

The temperature corresponding to the input pressure (vapor pressure) is found in TVAP. 

In addition, the enthalpy of the saturated liquid and the enthalpy of the saturated vapor are 
available through the parameter list; thus 

HGAS = HG, and 
HLIQ = HL . 

Tne sequence with the other saturation routine is as follows: 

T = y. y 

PVAP = TSATH(T, HG, HL) or PVAP = ETSATH(T, HG, HL). 

Here the vapor pressure is returned in PVAP. and again the enthalpies of saturated 
liquid anc vapor are available through the parameter list. Next, assume the call above 
has been made, but what really is required is the entropy of the saturated vapor. The 
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sequence is completed by 

SGAS - PHENTR (PVAP, HG) or SGASE = EHENTR (PVAP, HG). 

4. 3 Exceeding the Range of the Programs 
It is, of course, easy to violate the range of the programs when setting the values 
of the input variables. The problem seems to occur frequently if unit conversions are 
involved. We have omitted error messages for cases of boundary violations because each 
installation handles printed error messages differently. 

For the PT (ET) and PH (EH) programs, the out of bounds variable is automatically 
set to a value acceptable to the program, and computations proceed. The value returned 
is, of course, in error. Possible boundary violations are shown in the phase diagrams, 
figures 3 and 4, by dots, and the arrows indicate what conditions are actually processed. 

A complete description is as follows: 

Pressure P out of bounds 

If the entry pressure is less than 1 psia, the value of pressure is set to 1. 

If the entry pressure is greater than 5000 psia, the returned value is extrapolated from 
within the defined range. 

Temperature T out of bounds 

If the entry temperature i3 less than the temperature appropriate to melting conditions, the 
temperature is set to the melting temperature at the same pressure. 

If the entry temperature is greater than 6000° R, the value of temperature is set to 6000° R. 
Enthalpy H out of bounds 

If the entry enthalpy is less than the enthalpy appropriate to melting conditions, the enthalpy 
is set to the melting enthalpy. 

If the entry enthalpy is greater than 25, 000 BTU/lb, the value of the enthalpy is set to 
25, 000 BTU/lb. 

Saturation routines 

If either program is entered with a value outside the pressure or temperature range of the 
vapor pressure curve, the variable is set to the nearest end point of the curve (see figure 3). 
4. 4 Conversion to Other Systems or Installations 
The programs as listed in Appendix B were written to avoid known incompatibilities. 
We have tested them on a CDC 3600, a CDC 3800, and an IBM 360. Other installations have 
run them on an 1108 Univac. Thus, conversion requirements should be minimal. In- 
compatibilities that we are aware of include 


DATA / 


vs. 

DATA (. . . 

) 

ABSF ( 


vs. 

ABS ( 

. . > 

EXPF ( 

. . . ) 

vs. 

EXP ( 

...) 

Another source of 

incompatibilities 

arises from the 

use of different key board 

s on certain 

types 

of card punches, 

i. e. , ( ) + 

= vs. % < k 


%' ■« 


The programmer should note that the data arrays are entered via DATA statements. 
The upper limit on the number of such statements varies from installation to installation. 
The least number that we are aware of has been used. Since this means that in many cases 
not all of the data can be loaded with the actual program deck, we had recourse to data 
loading subroutines which are listed in tobies 4 and 5. If a limitation is encountered, a 
possible solution is as follows (Goldberg, 1967): 

Problem: File length not sufficient to include all DATA statements of routine 


PTENTH. 


Solution: Two segments used, 1. BLOCK DATA, 2. Function PTENTH modified. 

Details: 1. BLOCK DATA (new for PTENTH function) COMMON / B LKPTH/ AA. .. 

AO, followed by the DATA statements for arrays AA through AO (total 
1183 values equivalent to H array). 

2. Main function PTENTH is revised to include all DATA arrays 
except H (1183). The H array is linked to the main function through 
named common: COMMON/BLKPTH/H(l 183). 

4. 5 Pitfalls, Common Errors in Applications 
The earlier version of these decks has been available for nearly five years. 

Several points, where users have encountered difficulties, are described below: 

Storage Limitation - Each program is designed to be used independently of all others. 

Specify only those decks actually required. Total storage required is specified in tables 
4 and 5. Storage requirements range from 795 to 2151. 

Program Range - Observe the limits of the programs 1 - 5000 psia, 24. 845 - 6000° R, 

-130 - 25, 000 BTU/lb. Valid answers will not be obtained if these bounds are violated. 
Lifluid-Vapor Two Phase Conditions - There is no built-in warning in any of the PH(EH) 
decks telling the user that he is in the liquid-vapor two phase regie This means the user 
should check the value of Q to confirm the conditions he is trying to process. 

Liquid -Vapor T wo Phase Returns - Returns from this region for the variables k, 7), w, C 

p 

Vt And L, are artificial, as defined in table 3. The programs will return numbers for 
these variables, but what usually is desired are not these returns but rather the property 
values for saturated liquid or saturated vapor. 

E* ror * * The two categories of errors, interpolation error and estimated source data errors, 
are discussed in detail in Section 6. The user has to be aware of the limitations inherent 
in linear interpolation. The detailed graphs of maximum interpolation errors given in 
Appendix C should alert the user to possible problems, and serve as an aid in solving those 
problems. 

4. 6 Debugging Aids 

In addition to the program listings given in Appendix B, we have found that a set of 
selected values for the various properties is indispensable. In Appendix A we list values 


Details: 
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obtained from the PT(ET) programs. Values for all properties for a representative set 
of temperatures are printed for isobars of 1, 15, 50, 150, 500, 1500, and 4500 psia. 

These skeleton tables are to be used to establish approximate values, to make an order of 
magnitude check, to check units, and to debug the running of the programs at other 
installations. 

5. STRUCTURAL DETAILS OF THE PROGRAMS 
As mentioned before, each program can be visualized as an array of data coupled with 
an interpolation scheme. Each property value entered into the data array corresponds to a 
particular pair of P-T or P-H. The spacing or increment for each entry variable is 
independent of the other and varies over the surface of the property. The total network of 
spacings is called the grid and a sequence of equal spacings defines a region. 

5. 1 Grid and Region Layout 

The grid for each variable consists of 1 2 to 25 different regions. The size of each 
region and its relation to neighboring regions is determined by the desired interpolation 
error and by the relative linearity of the particular property surface. A plot of constant 
property lines is used as < *id in layout. An example of grid-region layout for p-hydrogen 
is shown in figure 5 where lines of constant C (in BTU lb- °R ) are plotted for entry 
variables P and H. Lines of constant C and regions G through Q are shown. Within each 
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region increments of the entry variables, P-T or P-H, are fixed. The spacmgs used for 
the various regions were intended to yield a 1% interpolation error; however, subsequent 
checks show that we have not always achieved this goal. A data array developed according 
to this scheme will vary from variable to variable; it will minimize the number of points 
required to represent the surface of the variable to be interpolated. 

While in general the property values entered in the data arrays are exact, there are 
a few cases where these entries have been adjusted to yield better results for the inter- 
polation. One example is entries at low pressures, that is, some regions are extended to 
negative pressures and the entries at these negative pressures are bogus values. They do, 
however, insure a more precise interpolation down to the lower limit of pressure. Region 
B shown in the partial grid layout of figure 6 is a typical example. The lower pressure limit 
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Figure 6. Partial Grid Layout for Program PTENTH 
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for interpolation from this region is 80 psia and entries in the data array at -100 psia are 
bogus values (shown as O). A similar technique better described as a partial overlap of 
region boundaries is used in Region A, C, and D. For best accuracy values at 15 psia are 
interpolated from Region D, not from Regions A, B, or C. Region E is the only one 
illustrated in which bogus values are not used. The entire sequence illustrates how difficult 
it is to follow the large variation in enthalpy in the region of dissociation. 

Interpolation along the two phase boundaries also involves the use of bogus points 
that have no counter part on the surface of the property in question. For the P-T programs, 
the points are across the two phase boundary while for the P-H programs they have been 
added inside the two phase boundaries. In either case, these bogus points are used only as 
interpolation aids to get from the single phase region right up to the two phase boundary. 

5. 2 Speed 

The programs gain their operational speed from performing two dimensional linear 
interpolation in a fixed array of data. Actual speed of computation depends, of course, on 
the particular installation. At our installation, a CDC 3800, returns run from 30, 000 to 
60, 000 per minute. The data used for interpolation have been placed into compact DATA 
statements to provide rapid access. Access to any given region is just as fast as to any 
other region because of the way searching is arranged. The first branching divides the 
number of regions in half, the second in half again, and so on. Speed is reduced by about 
a factor of 2 if lookup of the phase boundaries is required for the point under consideration. 
The programs return a single answer with two variables (three for PHTEMP) going 
through the transfer vector. 

5. 3 Continuous Functions 

In the approach taken here, continuous functions cannot be obtained because the 
increments of the entry variables change when crossing region boundaries. The problem is 
illustrated by two simple examples shown in figure 7. In case I, points 1, 2, 3, 4, and 5 lie 
on the boundary of the lower region, while points 1, 3, and 5 are points on the boundary of 
the adjacent upper region. Linear interpolation along the boundary for each region is shown 
by the straight lines. If points 2 and 4 are adjusted to 2' and 4', interpolation in both regions 
will be the same and jumps in the value of the property returned will not occur across this 
boundary. In case II, points 1, 2, 4, and 5 are in the lower region while 1, 3, and 5 lie in 
the adjacent upper region. The values of 2 and 4 cannot be adjusted to avoid a jump in the 
property interpolated*. An actual example where boundaries cannot be made to match can be 
seen in the transition from 30 psia, Region D, to 30 psia, Region C, of figure 6. 



\ 


* We cannot adjust the value of point 3 as the interpolation error at 3 would be larger than 
the desired 1%. 
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If the points on the boundary are common to both regions, then the 
boundaries can be made continuous. 


Figure 7. Matching of Values at Boundaries 
5. 4 Modifications, Future Developments 

Changes are a way of life with computer programs. Improvements will be dictated 

by the user — they may arise from better source data. One possible modification is 
described below. 

Return of Values on the Two Phase Boundary - In many applications involving the PH (EH) 
programs, some of the returns are from the two phase region. In this case, the property 
values on the respective boundaries are required. Since the programs calculate the values 
at the boundaries for the regular return, a slight change in program structure will yield 
the desired answers. The changes are illustrated by the following example for entropy: 

FUNCTION PHENTR (PRES, ENTH) change to 

FUNCTION PHENTR (PRESS, ENTH, SGG, SLL) 
and add at 25 + 2 SGG = SGAS 

SLL = SLIQ 

In other words, the new logic is similar to that used in PHTEMP(P, H. Q) and is accomplish* 
at the expense of having two additional variables pass through the transfer vector. 

Fu t ur e Possibilities - It is interesting to speculate on improvements that might be accompli: 
ed in the future based on what we have learned about the linear interpolation scheme. An ad 
vanced scheme of interpolation would save storage, improve accuracy, and allow a large 

combination of variables. It would be based on the following elements: 

2 




1 . 


The development of a single program to include 

a) an interpolation package 

b) a data array of fixed size, about 350 points per property. 

The use of higher order or, perhaps spline interpolation in the interpolation 
package. 

The use of the distribution in grid entries obtained during the present error 
analysis (see for example tables 8 and 9) to create an optimum point distribution 
for the new array of data. 

4. A separate program which would call the generating functions (see Section 6. 8) 
and create the array of data desired. 


6. ESTIMATION OF ERRORS 

Terms such as accuracy, precision, uncertainty, and error are used in scientific and 
engineering papers with many shades of meaning. To avoid ambiguity, the terms used in 
this report are defined as follows: 

First, we define error as the quantitative difference which exists between the true 
value and the value which is measured, calculated, or otherwise determined. 

Second, uncertainty is defined as an estimate of the error that may exist between 
reported and true values. Since usually we do not know the true value, the estimate is by 
nature a variable quantity. The estimate of error or uncertainty will vary according to the 
information available and the manner in which the estimate is obtained. 

Statistical methods normally are used to arrive at estimates of error or uncertainties. 
However, statistical methods are valid only for random processes, and therefore uncertain- 
ties based on statistical methods are estimates based on the presence of random errors only. 
In contrast to statistics, the nature of thermophysical data is such that systematic errors 
usually are present, i.e. , errors which are not random. Not only are systematic errors 
present but often they are of much greater importance than the truly random errors. 

To be realistic in assigning uncertainties we have combined the estimates of both 
random and systematic errors, except w^ere the uncertainties are labeled "statistical 
uncertainties" in which case only random errors are given. To minimize the risk of a true 
value being outside of the interval specified by our assigned uncertainty, the assignment of 
uncertainties is large enough to give a 99% chance of including the correct value, i.e. , a 
3 0 level. The a is used here in the conventional sense of standard deviation. 

Differentiation is made between uncertainties in the data of the sources chosen and 
interpolation error, i.e. , the difference between the values returned by linear interpolation 
and the source chosen. Maximum interpolation errors and estimated source errors for each 
program are displayed graphically in Append Jc C in three ranges of temperature: 25 to 
180 R, 180 to 1 200° R , and 1200 to 6000° R. The ranges of enthalpy for the PH and EH type 
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programs correspond closely to the temperature ranges just mentioned. All plots are given 
in per c entages exc ept for the low temperatur e range of enthalpy , i.e., 25 to 180°R for 
program PTENTH (ETENTH). The enthalpy in this range of temperature can assume both 
negative and positive values (see also fig. 4); therefore, the interpolation and source errors 

for this range only are plotted in BTU/lb. The s-rnbol in the plots signifies that the 

PVT surface for these temperatures is extrapolated. 

Estimated source errors are discussed in some detail in Section 7. The estimated 
source errors are established by the authors of the source data, by comparison with experi- 
mental values, by noting the spread of values if several experimental sets exist, or by 
interci mparison of several methods of calculation if no experimental values are available. 

For the high temperature extrapolations we combine the estimates of Woolley (1972) for the 
ideal gas, our estimates of errors for the real gas contribution, and an arbitrary estimate 

for errors due to dissociation. 1 

Interpolation errors are apt to intrude on most applications. These errors are 
nominally 3% except for the critical region and the region of dissociation where departures 

from be :t values become large. The critical parameters are P = 187.506 psia « 

C r ' j 

T = 59. 3 5 7° R , and H = 16.455 BTU /lb for p- hydrogen or H = 31.225 BTU /lb for 

C c 4 

e-hydrogen. The onset of dissociation varies from 2700°R at the lower pressures to 5400°R V 

at the highest pressures. The interpolation errors plotted are the envelope of maximum 
> error in interpolation at the temperature or enthalpy in question. They are established from 

an intercomparison of values returned by the interpolation routine and the generating 
functions. The comparisons are made at a large number of temperatures and a large number 
of pressures for each temperature, selecting the largest deviation in a given range of tem- 
perature or enthalpy for the plot. The distribution of temperatures and pressures at which 
the intercomparisons were made are established from the grids of the properties by 

collecting values at 1/2 of each normal grid increment and eliminating duplicates. As an 9 

example, each program was checked in the low temperature range at every combination of 
P and T given in tables 6 and 7. 


Table 6. Low Temperature Error Grid, Pressures, psia 


1,022 

36 . 74 J 

146.960 

225.000 

308.616 

41*. *36 

734.600 

1763.520 

1.500 

<♦ 4.088 

154.306 

230.000 

315.964 

426 . 1*4 

771. 540 

2000.000 

2.000 

51.436 

161. 656 

235.000 

323.312 

433.532 

806.280 

2500.000 

3.000 

58.784 

169.004 

240.000 

330.660 

440.680 

845.020 

2645.260 

4*000 

66.132 

176.352 

245.000 

338.006 

459.250 

661.760 

3233. 120 

5.000 

73.480 

160.026 

250.000 

345. 356 

477.620 

916.500 

3820.960 

6.000 

60.628 

183.700 

255.000 

352.704 

495.990 

955.240 

4406*600 

7*000 

86.176 

167.370 

260.000 

360.052 

514*360 

991.960 

5000.000 

a. ooo 

95.524 

190.000 

265.000 

367.400 

532. 730 

1026.720 


9.000 

102.872 

! 95. 000 

270.000 

374.746 

551.100 

1065*460 


10.000 

110.220 

200.000 

275.000 

362.096 

569.470 

1102.200 


11.022 

117.666 

205.000 

260.000 

389.444 

567.640 

1136.940 


14.696 

124.916 

210.000 

265.000 

396.792 

624.560 

1175.660 


22.044 

132.264 

215.000 

293.920 

404.140 

661.320 

1322.640 


29.302 

139.612 

220.300 

301.266 

411.466 

698.060 

1469.600 
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Table 7. Low Temperature Error Grid, Temperatures, °R 


25.000 

40.700 

55. 100 

26.000 

41.600 

56 . 00 G 

27.500 

42.500 

56.500 

29.000 

43.400 

57.000 

30.500 

44.300 

57.500 

30-900 

45.200 

58.000 

31.700 

46.100 

58.500 

32 . 6 U 0 

47.000 

59.000 

33.500 

47 . 9 UQ 

59.400 

36.400 

48.800 

59.625 

35.300 

49.700 

59.850 

36.200 

50.600 

60.000 

37.100 

51.500 

60.500 

36.000 

52. *400 

61.000 

36.900 

53.300 

61.500 

39.600 

54.200 

62.000 


6?. 500 

70. 50 0 

90.000 

63.000 

70.571 

S 3 • 600 

63.500 

71.500 

97.200 

64.000 

72 . COO 

100.800 

64.500 

72.900 

104.400 

65.000 

73.600 

108.000 

65.500 

74.700 

112.500 

66.000 

75.600 

117.000 

66.500 

76.500 

121.500 

67.000 

77.400 

126.000 

67.500 

78.300 

135.000 

68.000 

79.200 

144.000 

68.500 

60.000 

153.000 

69.000 

62.500 

162.000 

69.500 

85.000 

171.000 

70.000 

66.400 

180.000 


For the most part, users will be concerned only with the interpolation error because 
in most cases the uncertainties in the source data are negligible in comparison to the inter- 
polation error. Interpolation errors can be established quickly by visual estimate from the 
appropriate plot. If the generating functions (Section 7.8) are used, only the estimated 
source errors need to be considered. 


7. SOURCES OF PROPERTIES DATA 

To calculate values of the various properties, we require the thermodynamic functions 
of the ideal gas, a description of the PVT surface, a model .o calculate the effects of 
dissociation, and descriptions of the behavior of viscosity and thermal conductivity. The 
major references used are shown m figure 8 arid are discussed in more detail below. 

7. 1 Ideal Gas Properties 

As sources for the thermodynamic properties of the ideal gas we have considered 
Woolley, et al. (1948),NBS Circular 564 (Hilsenrath, et al., 1955) and NBS Monograph 20 
(Haar, et al., 1961). The first two sources are virtually identical; they differ only because 
Woolley includes the effect of nuclear spin in his tabulations for entropy (i.e. 2 R In 2). 

Normally the nuclear spin contributions are omitted; both Circular 564 and Monograph 20 are 
tabulated according to this convention. Slight differences in thermal properties between 
Circular 564 and Monograph 20 exist. They can be djtected at temperatures above 1260'R 
(700 K); they increase as the temperature increases. For C° the difference is about 1% at 
5000' R (2800 K) while for H° - Eg the difference is 44.6 BTU/lb (209 J/mol). A more 
recent calculation (Baehr, et al., 1966) is very close to the earlier values of Woolley, 
e’. al. (1948). We have, tnerefore, elected to use the properties given by Woolley for para- 
hydrogen up to 3600'R (2000 K). At this temperature Woolley discontinues the entries for 
parahydrogen but continues the entries for normal hydrogen on up co 9000' R (5000 K). At 
3600' R, differences between the ideal gas properties of para and normal hydrogen are 
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Figure 8. Primary Sources of Data 
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I negligible except in the case of entropy where the difference is 1.367 BTU/lb- # R. The 

j difference represents the change in order or configuration between para and normal hydrogen. 

Above 3600°R the next temperature entry in Woolley's tables is 5400* R (3000 K). For tern- 
peratures of 5400* R and above, the ideal gas properties of normal hydrogen were used for 
both para and equilibrium hydrogen since at these temperatures, equilibrium and normal 
hydrogen properties are identical. 

While the entries for equilibrium hydrogen are perfectly proper, for parahydrogen 
this means that we assume complete conversion to normal hydrogen between 3600* and 
5400 # R. A higher order interpolation insures that the transition in entropy for the para- 
j hydrogen decks is smooth and continuous. It becomes necessary to switch to the normal 

j mixture at some point because the dissociation calculations (Section 7.5) are based on 

normal hydrogen. It should be noted that the present programs are oriented toward a partic- 
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ular process which begins in the cryogenic range and ends at high temperatures. It would 
not be proper to use the parahydrogen programs on a process which begins in the cryogenic 
range, goes to high temperatures, and then is returned to the cryogenic range, because 
most likely the hydrogen would not be returned to the para form on cooling. 

7.2 The PVT Surface at Low Temperatures 

Values used in the data arrays for temperatures up to 180°R (100 K), for all prop- 
erties except viscosity and thermal conductivity, were obtained from NBS Monograph 94 
(Roder, et al., 1965). The uncertainties in source data were taken from this reference and 
are displayed as estimated source errors in Appendix C for temperatures from 25 to 180°R 
(-130 to 425 BTU/lb or 675 BTU/lb). 

7.3 The PVT Surface at Intermediate Temperatures 

Initially, we considered using the exponential equation of state (Woolley, et al., 1948). 
The difficulties with this equation are that it cannot be integrated in closed form to find the 
derivjd properties, it has a limited range, 500 to 1200° R (273 to 673 K), and the uncertain- 
ties in the derived properties cannot be calculated easily. 

As an alternative we considered a modified Benedict- Webb-Rubin (1940) equation of 
state. The advantages of an equation of state of this type are that the functions are internally 
consistent and continuous. In addition, statistical uncertainties for both PVT and derived 
properties can be obtained from the least squares estimation of the parameters. Various 
modifications of this equation have been used to correlate the PVT surfaces of nitrogen 
(Strobridge, 1962), helium (Mann, 1962), carbon monoxide (Hust and Stewart, 1963), neon 
(McCarty and Stewart, 1965), argon (Gosman, et al., 1969), and oxygen (Stewart, 1966). The 
equation has also been used earlier on parahydrogen (Roder and Goodwin, 1961) and on an 
inter comparison of normal and parahydrogen (Hust and Stewart, 1965). 

The PVT data used to estimate the parameters of the equation of state included the 
parahydrogen isotherms from 50 to 100 K (Gcodwin, et al., 1963), all of the PVT sources 
surveyed by Woolley, et al. (1946), and the isotherms of normal hydrogen reported by 
Michels, et al. (1959). The equation is as follows: 


P r P RT 4 p a (N t T 4 N. 4 N 3 /T 4 N 4 /T a 4 N B /T 4 ) 

+ p 3 (N, T a + N,T + N„) + p 3 (N, /T a + N 10 /T, + N„ /T 4 ) e NlT e * 
+ N lt p 4 T + p* <N 1S /T a + N 14 /T 3 + N,,/T 4 ) e Nl,P * 

+ N 1# p* . 

Th« coefficient, ere given in table 6 for P in atm, p in mol/i, T in K, and R = 0. 082053 5 
f-atm/mol-K. 
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Table 8. Coefficients of the Equation of State 


Ni 

= 

1. 1389049685 

X 

10- 3 

N 19 = 

-3. 7301781349 

X 

10 3 

N, 

= 

1.8072093722 

X 

10" 1 

N„ = 

1.0789473341 

X 

10 6 

n 3 

= 

-5. 309,7164419 

X 

10 1 

N 1S = 

8.0475741674 

X 

io- 

n 4 

= 

2.3690885344 

X 

10 3 

& 

(ji 

11 

1.6581404268 

X 

io- 

N b 

= 

-1.6306395805 

X 

10® 

It 

£ 

-1.8158230417 

X 

10 1 

n 9 

= 

3.2453695439 

X 

10 -8 

II 

1C 

4. 9623738040 

X 

10 2 

N , 

= 

-6. 3275640592 

X 

10"’ 

N ie = 

5.262288563 

X 

io- 1 

N, 

= 

2. 3115255859 

X 

10“ 3 

n 17 = 

-.0018 



N, 

= 

3.4448764044 

X 

10 1 






The equation is valid over the range of the data and was used from 180 to 1200“R. 
Unfortunately the major sources disagree considerably at the common boundary of 180‘R 
(see, for example, Roder. et al., 1965) and the resulting jumps in the properties are seen 
clearly in the plots of maximum interpolation error of Appendix C. A comparison of 
calculated statistical uncertainties for PVT and derived properties with the results of Roder. 
etal. (1965), Woolley, et al. (1948), and Michels, et al. (1959, 1963) is shown in figure 9 
where the base line is the present equation of state. The comparison confirms that the values 

obtained for a are realistic for both PVT and derived properties. For temperatures from 

180 to 1200-R and enthalpies from 425 to 5000 BTU/lb the calculated uncertainties are shown 
as source errors in the plots of Appendix C. 

7.4. The PVT Extrapolation 

The highest temperature for experimental PVT values is 673 K. This means the PVT 
aurface from 1200 to 6000»R has to be extrapolated. At these temperatures the gas is 
nearly ideal, that is, only the lower virial coefficients contribute to the equation of state. 

As a first step one has to know how the virial coefficient B varies as a function of temper- 
ature. and then how much the B(T) contribute, to the equation of state. The B for hydrogen 
attains it. maximum value near 300 K. Based on what is known for helium (Keesom, 1942). 
the B should decrease at higher temperatures. Woolley, et al. (1948) based the temperature 
dependence of B on this behavior and predicted the actual temperature dependence from the 
Lennard-Jone. (12-6) intermolecular potential function. Considering an extreme value, 
that is. 5000 p.i. (340 atm) and 5000-R (2800 K), the exponential equation of state (Woolley 
•t al.. 1948) yields a density of 0.337 lb/ft 3 (60 Amagat). In the range in which this qqu.tiol 
is valid. 360 to 1200-R ,200 to 673 K). the virial B contribute. 1... than 0.2% to the density 
for densities less than 60 Amagat. The contribution at higher temperature, can only be less; 
therefore, an extrapolation based on an equation of state which includes only the second 
virial coefficient is acceptable in terms of error. 
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Next we estimate the apparent error in the extrapolation of B(T). For 5400°R 
(3000 K), extrapolation of the Woolley, et al. (1948) calculation yields 14.55 cm 3 /mol for B 
while Fisher (1965) obtained 12.93 cm /mol. Fisher’s calculation is based on a '’softer” 
inter molecular potential function (i. e. , a 9-6) and he determines the molecular parameters 
from the high temperature viscosity measurements of Guevara and Wageman (1965). If 
hydrogen at this temperature behaved as an ideal gas, B would be zero; thus a reasonable 
estimate for the error in B at 5400°R (3000 K0 is ± 3 cm^/mol. 

For the data arrays we extrapolate the equation of state (Section 7.3). Since the 
equation was fitted over the entire range of data up to 1200°R (673 K), only one additional 
constraint (Hust and McCarty, 1967) was required to yield acceptable results. The equation 
wa ‘ unstrained to yield 14.55 cmVol for B(T) at 540C°R (3000 K). Z is approximately 
1.04 at 5400° R and 5000 psia, thus the real gas contribution is 0.04. Woolley (1972) 
estimates that the real gas contribution is uncertain by no more than 10%; therefore, we 
estimate the uncertainty in density to be 0.8% for these conditions. Uncertainties in the 
derived properties are estimated in a similar fashion assuming an uncertainty of no more 
than 10% for the real gas contribution, i. e. , the contribution from non ideality. This 
estimate assumes virtually no error in the ideal gas thermodynamic functions and virtually 
no error in the dissociation calculation (Woolley, 1972). 

7. 5 Dissociation 

At the highest temperatures considered in this report, molecular hydrogen will 
dissociate into atoms. Dissociation is both temperature and pressure dependent. To include 
the effects of dissociation several assumptions are required. 

!• Both ortho-para and dissociation equilibria occur instantaneously. 

2. Dissociation occurs according to the equilibrium constant for normal hydrogen. 

3. The property values can be obtained by adding the dissociation effects calculated for the 
id*al gas to the extrapolated properties of the real gas. 

Fortunately, the rate of ortho, para conversion is quite high (Farkas, 1935). It is 
proportional to the amount of atomic hydrogen. The dissociation rates are apparently also 
quite high (Sutton, 1962). Aseumption 1 is approximated, and thus assumption 2 is plausible. 
Hence, it is not necessary to consider the "hypothetical" dissociation constant for para- 
hydrogen (Farkas, 1935) even though the hydrogen is initially para. 

The dissociation constant is calculated from the ideal gas free energies. As indicated 
previously, the inclusion or exclusion of nuclear spin will affect entropy as well as free 
energy. The tables of atomic hydrogen from NBS Circular 564 (Hilsenrath, et al. , 1955) 
and the tables for molecular normal hydrogen from NBS Monograph 20 (Haar, et al. , 1961) 
were used. These tables according to convention.do not include the effect of nuclear spin. 
Equation 3.5 for free energy given by Woolley, et al., ( 1948 ) is equivalent to the table for 
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atomic hydrogen in Circular 564 but includes nuclear spin. The equations used to calculate 
dissociation effects of the ideal gas are those given by Woolley (1955). 

The effects of dissociation on a variable such as enthalpy are large. They are 
illustrated for a limited range of temperatures and a pressure of 50 psia in figure 10 (see 
also figures 11 and 12). Enthalpy differences are plotted vertically and the dissociation 
contribution ranges from 2% to 10% of the total enthalpy. Also shown is the linear inter- 
polation provided by program PTENTH(P,T). 



Figure 10. The Effect of Dissociation on Enthalpy 
7.6 Viscosity 

The primary experimental sources to consider are Diller (1965) for low temperatures, 
and Trautz and Zink (1930) and Guevara and Wageman (1965) for high temperatures. 
Differences in the latter two sources have now been resolved in favor of the more recent 
data (Guevara, et al. , 1971) allowing a critical evaluation of viscosity and thermal 
conductivity coefficients for the dilute gas (Hanley, et al. , 1970). The coefficients are 
represented by standard kinetic theory expressions and are based on a model inter molecular 
potential function. 
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For temperatures up to 180 °R for both dilute and dense gas, the analytical represen- 
tation of the results by Diller (1965) was chosen as source data. At intermediate and high 
temperatures including dissociation, the evaluation of Hanley, et al. , (1970) was chosen for 
the dilute gas while the dense gas contributions are calculated from the theory of Enskog 
(Hanley, et ai. , 1972). The dense gas contribution vanishes at 1000 K, that is. above 
1000 K the viscosity is independent of pressure or density but depends only on temperature. 
At yet higher temperatures with the onset of dissociation the viscosity becomes again 
pressure dependent as shown in figure 11 taken from Hanley, et al. (1970). We note that 

at the extreme temperature of 6000 ”R (3333. K) the pressure dependence is still relatively 
small. 

The uncertainty in viscosity is estimated as follows: About 2% for temperature up to 
1080'R, and then gradually increasing to about 5% at 6000°R. 
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Figure 11. The Viacoaity of Dissociating Hydrogen at Varioua Preaaurea 


7. 7 Thermal Conductivity 

For temperatures up to 180 *R,the analytical representation of the experimental 
results by Roder and Diller (1970) was chosen as source data for para hydrogen for dilute 
gas, dense gas, and liquid states. Dilute gas values for normal and equilibrium hydrogen 
differ considerably from the para values and are calculated from the change in the ideal 
gs.s heat capacity, C° v . The density dependence for n- and e- hydrogen is taken to be the 
same as para. 

At higher temperatures we again select the critical evaluation by Hanley, et al., (1970) 
for the dilute gas and the Enskog theory calculation (Hanley, et al. , 1972) for the density 
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(pressure) dependence. At 180 ° Retaking the difference between nornnal and para hydrogen 
into account, the agreement between the experimental results by Roder and Diller (1970, 
para), Golubev and Kalsina (1964, normal), and the calculation by Hanley, et al. (1972) is 
as good as 1% for pressures up to 1500 psia and no worse than 4% at pressures of 5000 
psia. Similar to viscosity, the dense gas contribution for thermal conductivity vanishes at a 
particular temperature (1000 K). At higher temperatures the potential function calculations 
agree reasonably well with the data of Blais and Mann (1960). At yet higher temperatures, 
the effects of dissociation on thermal conductivity become marked, see figure 12. The 



calculated, 

Hanley, et al. (1970) 

• • • 

experimental, 
Bethringer, et al. (1968) 


Figure 12. The Thermal Conductivity of Dissociating Hydrogen at Various Pressures 


calculations become sensitive to the exact model or intermoleculer potential chosen, see 
for example, figure 9 in the paper of Grier (1962). Two sets of data exist for the thermal 
conductivity in the region of dissociation. The theoretical estimates selected (Hanley, et al., 
1970) are in excellent agreement with the shock wave experiment of Bethringer, et al. (1968). 
The data of Israel, et al., (1966) differ markedly from the former experiment and from any 
potential function model that could be chosen. A possible explanation is that the sintered 
material used in this experiment enhanced dissociation. It almost appears as if the 
dissociation temperatures are lowered approximately 250 K by an apparatus effect. 
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The uncertainty in thermal conductivity is estimated to be 3% for temperatures 
between 25 and 180 ° R, increasing to 6% near 370 °R ; about 6% for temperatures up to 
3300 °R and then increasing gradually to 10% at 6000 °R. Note, however, that thermal 
conductivity is expected to be infinite at the critical point. 

7. 8 Generating Functions 

The programs described in the bulk of this report are designed to perform linear 
interpolations in arrays of data of the various properties. Each property value used in one 
of the arrays has to be computed in a separate operation and entered into the array. We have 
used several major computer programs to accomplish this task; to distinguish them from 
the linear interpolation programs, we call the major programs generating functions. An 
example of one such generating function is the program VALUES or the subroutine THERMO 
which was developed in the preparation of the tables presented in NBS Monograph 94 (Roder, 
et al. , 1965). Since a number of users have requested copie* of the generating functions, 
we list them below. 

Table 9. Generating Functions * 


Temperature Range 

Program Name 

Input 

Output 


25 - 180 'R 

program VALUES or 
subroutine THERMO 

PT 

P, S, H, U, Cp, Cy, 

(■&),- (g) T -- 



function CONC 

T r P » Cp 

k 



function VISCOS 

P, T 

11 

1 


program STAT 

T 

o - p composition and 
ideal gas functions, 
thermodynamic 

•V 

.1 

180*R and up 

program GO 

P, T 

p. S , H ,U , Cp. c v , 
Y, aj, k, ti 


* A recent revision incorporates the low temperature decks except STAT, an improved 
equation of state for high temperatures, and a new extrapolation to 10,000 psia into 
a single program HYPROP, see McCarty and Weber (1972). 
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Tables of Values for Selected Isobars 


Parahydrogen (Output from PT programs) .... 
Equilibrium Hydrogen 'Output from ET programs) 
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APPENDIX U. PROGRAM LISTINGS: CDC AND IBM 
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IBM 
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PTDENS 


PTENTR 
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PTENTR 


PTCOND 


65 

PTCOND 


PTVISC 


68 

PTVISC 


PTSOUN 


72 

PTSOUN 


PTCP, PTCV, 

PTGAMM 

75 

PTCP, PTCV, 

PTGAMM 

PTLFAC 
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PTLFAC 


TSATH 


83 

TSATH 


PSATH 


84 

PSATH 


PHTEMP 
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PHTEMP 


PHDENS 
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PHDENS 


PHENTR 
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PHENTR 


PHCOND 
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PHCOND 


PHVISC 
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PHVISC 


PHSOUN 
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PHSOUN 
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PHGAMM 
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PHLFAC 


ETENTH 
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ETENTH 
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ETDENS 


ETENTR 
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ETENTR 
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ETCOND 


ETVISC 
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ETVISC 


ETSOUN 
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ETCP, ETCV, 

ETGAMM 

137 

ETCP, ETCV, 

ETGAMM 

ETLFAC 


142 

ETLFAC 


ETSATH 


146 

ETSATH 


EPSATH 


146 

EPSATH 


EHTEMP 


147 

EHTEMP 


EHDENS 


150 

EHDENS 


EHENTR 


153 

EHENTR 


EHCOND 


157 

EHCOND 


EHVISC 


160 

EHVISC 


EHSOUN 


165 

EHSOUN 


EHCP, EHCV, 

EHCAMM 
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EHLFAC 
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3M»37«»%S6t4«0*S10f528) 1182 185. , 1 66 ,5 , 167 .25, 187 . 4687 5 , 157 . 506, 187 ,6385/ 

DA?A(JP«%»%»3»3t$r3»3»3f3'%r3f6r<»»tf*»5f5rllf6t3»3> DATA TS/24. 845,27 . 07, 29.81,33.07,36.18, 39.96 r<r<r. 12, 48. 33,51.97 ,54 

DAT A (MX* 2, 2 ,1,1, 3, 1,1, 1,1, 2, 1,4, 2,4, 2, 3, 3, 9, 4, 1,1) 179,56.72,57.80,58.57,58.99,59.18,59.29,59.34,59.353,59.356,59.4/ 

OAT A (IP* I, 8,-1., 0,1., -18., 0,1469.6, 1469.6, 1469.6, 1469. DATA LOC/1 , 21 ,37 , 46 ,55, 90 , 112, 1 34, 155 , 176 ,204, 2 16 , 276, 296, 320 , 356 
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IF ( IP .GT *HX ( N) > IP= MX l N) GO TO 33 

F=F P~ IP 7 DO 8 1 = 2,19 

FP=1.0-F IF ( P-PS ( I ) ) 10,9,8 

FT= (T-BT (N) )/0T (N) 8 CONTINUE 
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FUNCTION ETVISC (PRES, TEMP) FUNCTION ET VI SC (PRES , TEMP ) 

COMMON/E VI SCD/V ($36) DIMENSION P S ( 20 ) , TS (2 0 ) , J P ( 2 1) , MX ( ?1 ) , L OC ( 21 ) , BP ( 2 1 ) , DP ( 2 1 ) , &T ( 2 1 ) 

DIMENSION PS ( 20 ) , TS (28) , JP (21) » MX ( 21 ) , LOC (21) , BP( 21) , DP (21 ) , BT (21) 1 , DT ( 21) , V (536) 
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APPENDIX C. PLOTS OF MAXIMUM INTERPOLATION 
ERROR AND ESTIMATED SOURCE ERROR 
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EHLFAC 217 


\ 


177 


jTE'PPOLATJOfi ERROR SCARCE ERROR INTERPOLAT I ON ERROR SOURCE ERROR If JTERPOLAT I ON ERROR 



8 









J»CE ERROR INTERPOLATION EPPOR s:J?CE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR 






















ERRCR INTERPOLATION ERROR SCARCE ERROR INTERPOLATION ERROR sc j? C £ £ RR0 R It jTEPPOLAT I or j ERROR 






VRCE ERROR I f JTERPOL AT I Of J ERROR gC'v.'RCE ERROR INTERPOLATION ERROR SCARCE ERROR INTERPOLAT I ON ERROR 











-RCE ERROR interpolation error SOURCE ERROR interpolation error sciPCE ERROR interpolation error 



9 







.•RCE E**OR INTERPOL AT I Of J ERROR SCARCE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR 








: ERROR INTERPOLATION ERROR SCARCE ERROR INTERPOLATION ERROR 




-RCE ERROR INTERPOLATION ERROR SCARCE ERROR INTERPOLATION ERROR SOURCE ERROR Interpol at ion ERROR 


«■ 





ENTHALPY. BTu/lB. 


189 




























SC'jRCE ERROR INTERPOLATION ERROR 







8 





AT I ON ERROR SOURCE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR 



99 
























SOURCE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR 







O'JRCE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR 













SOURCE ERROR inteppol at i on error source ERROR interpolation error soupce ERROR interpolation error 


ETLFAC 





TEMPERATURE, RANKINE 


/ 


x 



207 



SOURCE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR 




















CE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR SOURCE ERROR INTERPOLATION ERROR 


1 




Et;Trt*i.Pf.B?u/uB. 


213 











ERROR INTERPOLATION ERROR SCARCE ERRCR INTERPOLATION ERROR SCARCE ERROR 








